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Abstract 

We study, both analytically and numerically, the cascade of failures in two coupled network 
systems A and B, where multiple support-dependent relations are randomly built between nodes 
of networks A and B. In our model we assume that each node in one network can function only 
if it has at least a single support node in the other network. If both networks A and B are 
Erdos-Renyi networks, A and B, with (i) sizes N A and N B , (ii) average degrees a and b, and 

> : 

(iii) c AB N B support links from network A to B and c BA N B support links from network B to A, 
we find that under random attack with removal of fractions (1 — R A )N A and (1 — R B )N B nodes 
respectively, the percolating giant components of both networks at the end of the cascading failures, 
/x^ and jU^,, are given by the percolation laws = R A [1 — exp (— Cq — exp (— afj,^)] and 

^ = R B [1 - exp (-Cq B ^)][1 - exp (-bfi^)]. In the limit of c BA -> oo and c AB -> oo, both 

•rH . 

networks become independent, and the giant components are equivalent to a random attack on a 
single Erdos-Renyi network. We also test our theory on two coupled scale-free networks, and find 
good agreement with the simulations. 
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I. INTRODUCTION 



In recent years, there has been extensive effort to study and understand the properties of 



'ocused on properties of single networks which do 



complex networks. Research has mainly 
not interact or depend on other networks l Ml4|. Recently, the robustness of two interdepen- 
dent coupled networks has been studied [l5, [l6]. In interdependent networks, the failures 
of nodes in one network, A, will cause failures of dependent nodes in the other network, B, 
and vice versa. This process occurs recursively, and leads to a cascade of failures. It has 
been shown both analytically and numerically that the robustness of two interdependent 
networks is significantly lower compared to that of a single network [15j. Furthermore, the 
percolation transition in coupled networks is first order compared to the known second order 
transition in a single network jl^, 

Previous studies on two interdependent coupled networks are restricted by the condition, 
that to function each node in network A depends on one and only one node in network B 
and vice versa 15(. However, in the real world, this assumption may not be valid. A single 



node in network A may depend on more than one node in network B and will function as 
long as one of the support nodes in network B is still connected. Similarly, a node in network 
B may depend on more than one support nodes in network A. As long as one of the support 
nodes functions, the node in network B will also function. 

Examples of such systems include the coupled power grid network and the communication 
network which controls the power grid, where both networks depend on each other. In 
general, one power station provides power to more than one communication stations, and one 
communication station controls more than one power stations. As long as a communication 
station can obtain power from one power station, it can function properly. On the other 
hand, one communication station is sufficient in sending control signals which make one 
power station function properly. However without any power, the communication station 
will fail, and without control the power station will stop working. Indeed in the 2003, due 
to failure of some power stations in Italy, the communication control system was damaged. 
This damage caused further fragmentation of the power grid, which finally led to a blackout 
in a sinificant part of Italy. 

Under random attack, which is characterized by random removal of nodes in one or both 
networks, the coupled network systems demonstrate significantly different behavior from 
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that of a single network 15]. The failures of nodes in network A can lead to the failures 
of dependent nodes in network B, and the failures of nodes in network B can produce a 
feedback on network A leading to further failures in network A. This process can occur 
recursively and can lead to a cascade of failures. 

We provide a theoretical framework for understanding the robustness of interdependent 
networks with a random number of support and dependent relationships. Our theory agrees 
well with the numerical simulations of several model network systems, including coupled 
Erdos-Renyi (ER) 17] and coupled scale-free (SF) {2] networks. Our work extends previous 
works on coupled networks [3, from one-to-one dependent-support relation to multiple 
Poissonian dependent-support relation. Our model can help to further understand real-life 
coupled network systems, where complex dependence-support relations may exist. 

We define the stable state to be the state when the cascade of failures ends. We show 
that for two coupled ER networks the giant components of both networks in the stable state 
follow a simple law, which is equivalent to random percolation of a single network in the 
limit of a large number of support links. Our theory is relevant to a broad class of real-world 
interdependent network systems. 

The paper is organized as follows. In Sec. II, we explain the model of the cascade of 
failures with random support- dependent relations. In Sec. Ill, we derive analytically the 
process of the failure cascade. In Sec. IV, we present numerical tests on coupled ER and 
SF networks. 



II. THE MODEL 

We assume two networks A and B of sizes N A and N B and with given degree distributions, 
P (k) and P B (k), of "intra-links" connecting nodes in the same network (Fig. 1). The 
dependency relation is represented by a link connecting the support node in one network and 
the dependent node in the other network ("inter-links"). The inter- links between network 
A and network B are random and uni-directional. Initially (stage n = of the failure 
cascade), there are c BA N A inter-links distributed randomly from nodes in network B to 
nodes in network A, representing the dependencies of node in network A on nodes in network 
B. Similarly, there are c AB N B inter-links from nodes in network A to nodes in network B, 
representing dependencies of nodes in network B on nodes in network A. The dependent- 
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support relations are random i.e., for each inter- link, the support node and dependent node 
in the two networks are chosen randomly. c BA and c A b are the initial mean degrees of the 
corresponding inter-links for networks A and B respectively. 

In our model, one node in either network A or network B can have zero, one or several 
support nodes in the other network. We assume that to function each node in network A 
requires at least one support node in network B, and vice versa. Nodes in each network 
which are not connected to the giant components of that network by the intra-links, and 
nodes without support inter-links, are considered to be not functioning and regarded as 
failed nodes. 

The attack on the network is represented by a random removal of fraction 1 — R A nodes in 
network A and 1 — R B in network B, where in general, R A 7^ R B '. The process of the cascade 
of failures is demonstrated in Fig. 1. With N A = N B = 7, we show the case of random 
removal of one node in network A (R A = 6/7) and one node in network B (R B = 6/7). At 
each stage of the cascade of failures, both networks will experience further failures. Without 
loss of generality, we assume the random attack on network A occurs before that on network 
B. Thus, when we analyze the first stage of the cascade of failures in network A, all the 
support inter-links from network B are considered functioning. At each stage, the nodes 
which do not have any support inter-links from the other network, and the nodes which 
are separated from the giant component of the network, are considered to have failed. This 
process will continue until no further node failure in either network occurs. 

III. ANALYTICAL SOLUTION 

The stable state of the two stable connected giant components in both networks are 
usually reached after several stages in the cascade of failures. We define c AB and c BA as the 
average number of support inter-links remaining in the stage n of the cascade of failures of 
network A and network B respectively. 

In each stage of the cascade of failures, we analyze first network A then network B. 
Such a procedure does not have any affect on the final result of the failure cascade. While 
considering network A at stage n, we assume all the inter- links {c BA x ) from network B to 
network A of the previous stage are working. When considering network B at stage n, we 
use the updated c AB after the nth stage of network A. 
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FIG. 1: (Color online) Demonstration of the stages of the cascade of failures in coupled network 
A and B of size N A = N B = 7. Curves represent intra-links within the network, while arrows 
(directed links) represent the inter-links connecting a support node in one network to the dependent 
node in the other network. Among the total 12 directed links, half of them c BA N A = 6 represent 
the support from nodes in network B to nodes in network A (arrow from nodes in network B to 
nodes in network A). The rest c AB N B = 6 links represent the support from nodes in network A to 
nodes in network B. The support-dependent relations between nodes in network A and network B 
are random. Initially, the attack is on node A\ (shown in red) in network A and node Bq (shown 
in red) in network B. The failed nodes are removed from the plot. In the first stage of the cascade 
of failures in network A, A\ fails because of removal, node A7 fails because of it has no support 
inter-links, A2 and Aq fail because of separation from the giant component of network A. All the 
failed nodes will lead to failures of inter-links starting from them. Since we assume that the attack 
on network A occurs before that on network B, the inter-link from Bq to A5 is considered to be 
functional. In the first stage of the cascade of failures in network B, we first remove inter-links 
connecting network B to non-giant-component nodes in network A (B3 to Ax, B2 to A2, B-j to 
A7). Next, node Bq fails because of the attack and nodes Bi, B<i and B7 also fail because of no 
support. B3 fails because it becomes separated from the giant component (nodes -B4 and B5) of 
network B. Finally, after removing inter-links connecting from network A to non-giant-component 
nodes in network B (.A3 to B3, A4 to B3 and A5 to Bq), the coupled network system reaches a 
stable state after one step in the cascade of failures, since all nodes in both giant components are 
connected and each node have at least one support node from the other network. 

5 



Since, in our model we randomly set up c AB N B inter-links from network A to network B 
and c BA N A inter-links from network B to network A, the degree distributions of inter-links 
of both networks follow Poisson statistics. There exist nodes in network A and network B 
without support initially (nodes A-j and B\ in Fig.l). These nodes are regarded as failed 
and need to be considered in addition to the random removal of nodes. 

At stage 1, after removal of a fraction 1 — R A of nodes in network A, taking into account 
also the nodes with zero inter-links at stage 0, the result is equivalent to a random removal 
of a fraction 1 — p A (R A ,c BA ) of nodes in network A. The giant component of network A 
will constitute a fraction g A (pf, c BA ) of the remaining p A N A nodes of network A. Thus, the 
fraction /i A of the giant component with respect to the N A , the original size of network A, 
is 



Since we assume that all support links from network B exist, we use c B . 

The failures of (1 — fi A )N A nodes in network A will lead to failures of support inter-links 
from network A to network B. The working support inter-links from network A to network 
B, c AB , will decrease from c AB to c AB (fx A ). 

Next, at the first stage of cascade of failures in network B, a fraction 1 — R B of nodes in 
network B will fail because of the initial attack. Combining this fraction with the fraction 
of nodes in network B having zero degree of support links after the first stage of cascade 
of failures in network A, the joint effect is equivalent to a random removal of a fraction 
1 — pf(R B , c AB ) of nodes in network B. The giant component of network B will constitute a 
fraction gf (pf ,c AB ) of pfN B nodes. Thus the fraction of the giant component of network 
B after the 1st stage of cascade of failures is 



The number of support links c BA then will be reduced from c BA to cf j4 (/if ). 

The 2nd stage of the cascade of failures in network A is equivalent to the first stage 
in network A with the updated cf A replacing c BA . Accordingly, the 2nd stage of cascade 
failures in network B is equivalent to the first stage in network B with the updated c AB . 
During the cascade of failures, c AB and c BA decrease as n increases. 

In general, the nth stage of the cascade in network A is equivalent to the 1st stage of the 
cascade in network A with c BA replaced by c BA v Similarly, the nth stage of the cascade in 



^ A = p A (R A ,c BA )g A (p A ,c BA ). 



(1) 



^=P?(R B ,cf B )g B (p B ,c AB ). 



(2) 
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network B is equivalent to the 1st stage of the cascade in network B with c AB replaced c AB . 
The general forms of the giant components of both networks at stage n of the cascade of 
failures can be expressed as 

.a 



u B 

H'n 



(3) 



Pn 

where p A and p B are the equivalent fractions of nodes in network A and network B respec- 
tively after random removal, and g A and g B are the fractions of the giant components in the 
remaining p A and p B fraction of nodes. The key to the analytical solution of this process is 
to find the way how c AB and c BA decrease with the cascade stage n. 

Next, we will use the apparatus of generating functions 
of fi A and (1%, c AB 

network A and P B (k) of network B are 



„ , and c„ . The generating functions of the degree distribution P (k) of 



20j to derive the analytical forms 

A i 



G A0 (x) = ^ k=Q j 
G B0 (x) = Z%L P B (k)x k . 

Analogously, the generating functions of the underlying branching processes are 



(4) 



(5) 



G A1 (x)=G' Ao (x)/G' AO (l), 
G B1 {x) = G' m {x)/G' m {l). 

After random removal of a fraction 1 — p of nodes, the remaining p fraction of the network 
will have different degree distribution. The new generation functions Go and G\ will be 

HQ 

G A0 (x,p) = G A0 (1 -p(l - x)), 
G B o(x,p) = G B o{l - p(l - x)), 
G A1 (x,p) = G A i(l -p(l - x)), 
G m (x,p) = Gbi(1 -p(l - x)). 
Randomly connecting c BA N A support links from network B to network A, the degree 
distribution of the inter-links in network A follows a Poisson distribution with average degree 



(6) 



^BA 



Similar, for network B, 



P\k) 
P B (k) 



[c 



BAlk 



k\ 



C 



ABlk 



k\ 



(7) 
(8) 



During the process of the cascade of failures, because the support-dependent relations are 
uncorrelated with the network properties of network A and network B, P A {k) and P B {k) 
will remain Poisson distributions with the new c BA and c AB , which decrease as n increases. 

Initially, there will be a fraction P A (0) = e~ cBA of nodes in network A which do not 
have any support links from network B. Since the attack on the 1 — R A fraction of nodes 
from network A is random, there will be overlap between the attack and the initially not 
working nodes (without support links) in network A. The joint effect is equivalent to a 
random removal of a fraction 1 — R A (1 — e~ cBA ) of nodes in network A. What happens in 
network B is similar to network A with R A replaced by R B and c BA replaced by c AB . Thus 

p A = R A (l-e-^), 
p B = R B (l-e-^ B ), 

where c B _ x and c AB are average degree of inter-links of network A and network B at the end 
of stage n — 1 and n respectively. 

According to the results on single networks jisl . \\\ . after random removal of a fraction 
1 — p A (or 1 — p B ) of nodes, the fractions of nodes that belong to the giant components of the 
remaining network A or network B, which have p A and p B fractions of nodes respectively, 
are 

' 9 A (P A ) = l-G A0 (f^p A ), 
g B ( P B ) = l-G B0 (f B ,p B ). 

where f A and f B satisfy transcendental equations 

f A = GAl(fn,Pn), 
_ fn = GBl(fniPn)- 

Thus \i A and fi B , the fractions relative to their original sizes of giant components of network 



A and network B 



18 



19| are 

= Pu9 A (P A 



n ) i 



(12) 

»n =Pn9 B (Pn)- 



Accordingly, c AB and c BA follow the relations 



rAB _ „ab . A 
r BA _ „BA. B 



When the cascade of failures process stops, f A , , p A , p„ , 



AS BA 
5 c n 5 



and \i B all reach 



the constant values, f A , f^, p^, p^, , c£, , and /z^ respectively. In principle, these 
final values can be found from the set of equations: 

Zoo — ^l(/oo)Poo)) 
f B — GBl{fx,P^o), 

-C BA \ 

e °° ), 



Pi 



R A (1 

R B (l-e- c ~), 



^AB 



(14) 



Moo 



pig A {p. 



A ) 



^=P^9 B (Px)- 



The functional forms of Gai, G\bi, G^o and G^o can be complicated, thus Eqs. (|T4l) can 
only be solved numerically for most cases, including coupled SF networks. However, for ER 
networks, Gq(x) and G\(x) have the same simple form 20], 



G (x) = G 1 (x) = e W x - 1 \ 



(15) 



where for network A, (k) = a and for network B, (k) = b. Thus, the above process of the 
cascade of failures can be significantly simplified. Eqs. ffTOj) can be reduced to 



(16) 



g A {p A ) = l-f A , 
9 B (P B ) = 1 - f B - 

Excluding p^, p^, c^ s , c^ 4 , /i^ and /i^ from Eqs.( fT4l) . for the stable state of two coupled 
ER networks, we get a system of two equations with two remaining unknowns, f A and 



fi = G A i{fLR A {l-e-^ A ))=e 



aR A (f£-l)(l-e~^ ) = aR A (f A -l)(l-e% 



3 W? B (/£-l)(l-e c C/oo-1)) 



(17) 



fB = G Bl {f^R B {l - e" c ~ )) = eW-OS-Dd-.— ) , 

The solutions of Eqs. fllT]) for f A and can be used to calculate the stable mutually 
connected giant components of both networks as 



/4 



pig A (pi 



R A (1 



l )(l-f A ) = -Hf A )/a, 



(18) 



J&7 B (P&) = R B a ~ e- c ~ )(1 - f B ) = ~Hf B )/b. 
From Eqs. ljTB"]) . using Eqs. (|T3|) . we can derive the percolation law for the fractions of the 
stable giant components of both coupled ER networks 
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fJ ,A=R A (l- e - c Z ^)(l-e- a ^), (19) 
fj,g = R B (l-e- c $ B ^)(l-e- b ^). (20) 

Eq. (119p and Eg. ( 120]) are simple and can be related to the theory of random percolation 



of a single ER network [21|, |22j, for which the fraction of the giant component is /i c 
R(l — e - ^ M °°). The coupled ER networks bring new terms 1 — e~ c ° Afi ™ and 1 — e~ c " B ^°°. 
In the limit of c AB — > oo (or c BA — > oo), the giant component of network B (or network A) 
does not depend on the other network and behaves similarly to the random percolation of a 
single network. 



From Eqs. (fT9]) and ( 1201) . we find and for a given set of parameters R A , R B , c 



BA 



and c AB . However, for some values of R A , R B , c BA and c AB , the solutions for /i^, and /i^ 
between and 1 may not exist. There exist critical thresholds of R A , R B , c BA and c AB above 
which the two coupled ER networks have non-zero mutually connected giant components 
(see FigfJ]). They are represented as R A , R B , c BA and c AB . These values can be solved by 
finding the tangent point of the two curves is plotted as a function of fi^ as shown in 
Fig [5]) represented by Eqs. (fT9i) and ( 1201) . The thresholds can be found from the tangential 
condition 

together with Eqs.® and flU}. 



1, (21) 



IV. NUMERICAL SIMULATIONS 



Next, we compare our theoretical results obtained in Sec. Ill to results of numerical 
simulations. We begin with comparing the simulations of the stages of the failure cascade 
in coupled ER networks with our theoretical predictions. In all our simulations, we use 
N A = N B = 10 6 . Fig. [3] shows fi A and fi B as a function of n for a = b = 4, c AB = c BA = 4, 
R B = 1 and for different values of R A . One sees very good agreement between the theory 
and the simulations. Close to R A , both \i A and fi B show large fluctuations between different 
realizations (shown in Figs. 3c and 3d). The random realizations split into two classes: 
one that converges to a non-zero giant component for both networks and the other that 
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FIG. 2: Demonstration of the functional relation between /j,^ and //^ in Eq. (|19p and Eg. (|20j) for 
a system of two coupled ER networks with a = b = 4, = c BA = 4, and i2 B = 1 with different 
values of R A . Since we use R B = 1, at different R A , the relation between /j,^ and //^ given by 
Eq. (j20p remains the same (shown by the dashed line). Eq. (jl9p with R A = 0.6, 0.43 and 0.4 are 
shown. One can see that when R A < 0.43, there exists no solution (/U^,/^) for Eq, (jl9p and 
Eq.J2QI|. The value 1 - R A = 1 - 0.43 = 0.57 represents the maximum fraction of nodes in network 
A one can randomly remove at the initial stage of the cascade of failures for which the non-zero 
giant components of both networks still exist at the stable state. The abrupt fragmentation of the 
stable giant components at R A < R A represents the first order nature of the percolation phase 
transition. 

results in a complete fragmentation. The agreement between the simulations and theoretical 
predictions is also good for different values of R B , a, b and c AB and c BA . 

In Fig. HJ we compare the theoretical predictions and simulations of the giant components 
at stages of the cascade of failures for a system of two coupled SF networks with X A = \ B = 
2.5, c AB = c BA = 4, R B = 1 and different values of R A . Similarly, we obtain agreement 
between the theoretical predictions and the simulations. Close to R A , both \i A and \i B of 
different realizations show large fluctuations and the random realizations also split into two 
classes. We also simulated other values of R B , A A , X B , c AB and c BA and found very good 
agreement between and theoretical predictions and simulations. 

The fractions of the giant components of both network A (/i A ) and network B (/x A ) m the 
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Eq.(20) 



Eq.(19)R =0.6 
-Eq.(19)R A =0.43 
-Eq.(19)R A =0.4 




FIG. 3: The case of coupled ER networks. Comparison between the theoretical predictions, ob- 
tained from Eqs. (|16p . and Eg s, (| 121113]) . and numerical simulations with N A = N B = 10 6 , a = b = 4, 
c AB = c BA = 4, R B = 1 and several values of R A . (a) and (b) show \i A and \i B at different stages n 
of the cascade of failures for R A = 0.7 and 0.6 above R A ~ 0.43 for both theory (lines) and simula- 
tions (symbols). One can see that both \i A and \i B approach a stable value and at the end 
of the cascade of failures. The agreement between theory and numerical simulations is very good, 
(c) and (d) show fi A and fi B at different stages n of the cascade of failures for R A ~ R A . The bare 
lines represent several realizations of the simulations and the lines with symbols represent the the- 
oretical predictions. One can see that for the early stages (small n) the agreement is good, however 
at large n the deviation due to random fluctuations in the actual fraction of the giant component 
starts to increase. The random realizations split into two classes: one that converges to a non-zero 
giant component for both networks and the other that results in a complete fragmentation. 
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FIG. 4: The case of coupled SF networks. Comparison between the theoretical predictions, obtained 
from Eqs.JTD}, and Eqs.(|llJ13J), and numerical simulations with N A = N B = 10 6 , \ A = \ B = 2.5, 
c AB = c BA = 4, R B = 1 and different values of R A . (a) and (b) show [i A and \i B at different 
stages n of the cascade of failures for R A = 0.7 and 0.6 above R A ~ 0.385 for both theory (lines) 
and simulations (symbols). Similar to Fig. [31 one can see that both /i A and \x B approach a stable 
value /i^, and /x^ at the end of cascade failures. The agreement between the theory and numerical 
simulations is very good, (c) and (d) show /j A and [i B at different stages n of the cascade of failures 
for R A ~ R A . Bare lines represent several realizations of the simulations and the lines with symbols 
represent the theoretical predictions. One can see that for the early stages the agreement is good, 
however at large n the deviation due to random fluctuations in the actual fraction of the giant 
component increase. The random realizations split into two classes: one that converges to a non- 
zero giant component for both networks and the other that results in a complete fragmentation. 
For coupled SF networks, at R A = R A , the fluctuations of both \x A and \x B seem to be relatively 
larger than that of coupled ER networks, due to the existence of large degree nodes in SF networks. 
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FIG. 5: (a) /x A and (b) in the stable states as a function of R A and i? B for coupled ER networks 
A and B with N A = N B = 10 6 , a = b = 4, c AB = c BA = 4 and 1 - i? A = 2(1 - R B ). Several curves 



for c, 



4, 8 and 32 are shown. The theory (lines) fits very well with the simulation results 



(symbols). One can see that for a given set of c AB and c BA there exist critical thresholds R A and 
R B , below which both networks will collapse and have no stable non-zero giant components. The 
value of R A approaches the critical threshold of random percolation (r = 1/a = 0.25) of a single 
network for large values of c AB and c BA . The initial attack on network B is smaller than that on 
network A and thus R B > R A . 

stable state can be found from Eqs.(|T9l) and ( 120|) for coupled ER networks. We solve these 
equations numerically for different R A and R B , and compare the theoretical predictions with 
the simulation results (Fig. [5]). For simplicity, we assume a = b = 4 and that the initial 
fraction of nodes affected by the random attack in network A is twice as large as that in 
network B (1 - R A = 2(1 - R B )). We test different average degrees of inter-links for both 
networks c AB = c BA . 

In Fig. [5l we present results for the giant components of both networks as a function of 
R A and R B . For different sets of c AB and c BA , we find again that the theory fits well with 
simulation results. One can see the critical R A and R B , which are the minimum fractions 
of both networks needed to be kept at the beginning of the cascade of failures in order to 
have non-zero connected giant components of both networks at the stable state. At R A and 
R B , both /x A and //^ show an abrupt change from a finite fraction (n A and /if) to zero. As 
c AB and c BA increase, R A approaches the critical threshold of random percolation of a single 
ER network, which is 1/a. As expected for single networks, \i A and \i B approach and a 
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FIG. 6: The dependences of R c (a) and \i c (B) on c A = c B for coupled ER networks with 
N A = N B = 10 6 , a = b = 4, and 1 - R A = 2(1 - R B ). In (a), the critical initial fraction of the 
network A, R A , and the critical initial fraction of the network B, R B , are shown as a function of 
c AB = c BA . The theory (full line and dashed lines) fits well the simulation results (symbols). R A 
approaches the critical threshold (1/a) of random percolation for a single network, as predicted by 
Eqs. (|19p and ()20p . In (b), the giant component of both network \x A and /i B are shown as a function 
of c AB — c BA at R A and R B . The theory (full line and dashed lines) fits well the simulation results 
(symbols). One can see that for large c AB = c BA , \i A and fi B both approach zero as expected 
for a single network. However, \x A and \i B will never reach zero for finite c AB and c BA , the phase 
transition thus remain a first order. 

second order phase transition exists for large c AB and c BA . However, for finite c AB and c BA 
the changes of /i^ and fi^ are not continuous at R A and R B , indicating a first order phase 
transition. This result is predicted by Eqs. (I19j) and f )20|) . We find that the theory fits well 
with the simulation results for the entire range of R A and R B for different values of c AB and 
c . 

Next, we study the dependence of R A and R B , \i A and /if on c AB and c BA (see Figj6]). 
For simplicity and for comparing with our earlier cases, we use the same set of parameter 
for both network: a = b = 4, c AB = c BA , and 1 — R A = 2(1 — R B ). As seen from Figj6j 
the theory fit well the numerical simulations. For large c AB = c BA , one can see that R A 
approaches the random percolation threshold 1/a on a single ER network. This behavior 
indicates that when one network A has enough support from network B and vice versa, both 
networks will behave as if they are independent. Indeed for large c AB and c BA , at R A and 
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R B , the stable giant components of both networks \i A and \i A approach zero as expected for 
a second order percolation phase transition. However, as seen from Eqs. (I19j) and (1201) . for 
finite values of c AB and c BA , neither \i A nor \i B is zero. This result supports the existence of 
a first order phase transition for the entire range of c AB and c BA . Good agreement between 
theory and simulations, and similar behavior of R A , R B , \i A and \i A as a functions of c AB 
and c BA have been found for other sets of parameters. 



V. CONCLUSIONS AND DISCUSSIONS 

In this paper, we extend previous works [15I . \v\ on the cascade of failures on interde- 
pendent networks by considering random support-dependent relations between two coupled 
network systems. Our theory is in excellent agreement with the numerical simulations on 
coupled Erdos-Renyi (ER) and coupled scale-free (SF) networks systems. For coupled ER 
networks, the percolation law for the giant components of both networks have a simple 
form, which in the limit of large number of supports gives the percolation law for single 
networks. Only in the limit of a large number of support is the percolation transition of 
second order, while in general case, the coupled network show a first order phase transition. 
Our model can help to further understand real-life coupled network systems, where com- 
plex dependence-support relations exists. Recently, a complementary approach to study the 
robustness of coupled networks system has been proposed [23|, which is based on a quite 
different assumption about the way networks are coupled. In contrast to our case where p c 
increases due to coupling, in their case p c decreases. Note that there are also recent efforts 
to study the robustness of single networks [24M26I] undergoing targeted percolation, which 
correlates with the topology of the network. In the same spirit, our work can be extended 
to study the robustness of coupled networks under non-random percolation. A first attempt 



in this direction for interdependent networks can be found in Ref. 



23]. 



Acknowledgments 

We wish to thank the ONR, DTRA, EU project Epiwork, and the Israel Science Foun- 
dation for financial support. S.V.B. thanks the Office of the Academic Affairs of Yeshiva 
University for funding the Yeshiva University high performance computer cluster and ac- 



16 



knowledges the partial support of this research through the Dr. Bernard W. Gamson Com- 
putational Science Center at Yeshiva College. 



17 



[1] D. J. Watts and S. H. Strogatz, Nature (London) 393, 440442 (1998). 
[2] A. L. Barabasi and R. Albert, Science 286, 509(1999). 

[3] R. Pastor-Satorras and A. Vespignani, Evolution and Structure of the Internet: A Statistical 

Physics Approach (Cambridge University Press, 2006). 
[4] S. N. Dorogovtsev and J. F. F. Mendes, Evolution of Networks: from Biological nets to the 

Internet and WWW (Oxford University Press, New York, 2003). 
[5] R. Cohen, K. Erez, D. ben-Avraham, and S. Havlin, Phys. Rev. Lett. 85, 4626 (2000). 
[6] L. A. Braunstein, Z. Wu, Y. Chen, S. V. Buldyrev, T. Kalisky, S. Sreenivasan, R. Cohen, E. 

Lopez, S. Havlin, and H. E. Stanley, International Journal of Bifurcation and Chaos 17, 2215 

(2007). 

[7] J. Shao, S. Havlin, and H. E. Stanley, Phys. Rev. Lett. 103, 018701 (2009). 
[8] M. E. J. Newman, S. H. Strogatz, D. J. Watts, Phys. Rev. E 64, 026118 (2001). 
[9] M. E. J. Newman, Networks: An Introduction, (Oxford University Press, USA 2010). 
[10] R. Cohen and S. Havlin, Complex Networks: Structure, Robustness and Function, (Cambridge 

University Press, Cambridge, 2010). 
[11] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, D.-U. Hwang, Physics Reports 424, 
175(2006). 

[12] G. Caldarelli and A. Vespignani, Large scale structure and dynamics of complex networks, 

(World Scientific, Singapore, 2007). 
[13] A. Barrat, M. Barthelemy, and A. Vespignani, Dynamical processes on complex networks 

(Cambridge University Press, England, 2008). 
[14] D. S. Callaway, M. E. J. Newman, S. H. Strogatz and D. J. Watts, Phys. Rev. Lett. 85, 

5468(2000). 

[15] S. V. Buldyrev, R. Parshani, G. Paul, H. E. Stanley, and S. Havlin, Nature (London) 464, 
1025(2010). 

[16] R. Parshani, S. V. Buldyrev, and S. Havlin, Phys. Rev. Lett. 105, 048701(2010). 
[17] P. Erdos and A. Renyi, Publ. Math. 6, 290 (1959); Publ. Math. Inst. Hung. Acad. Sci. 5, 17 
(1960). 

[18] J. Shao, S. V. Buldyrev, R. Cohen, M. Kitsak, S. Havlin, and H. E. Stanley, Europhys. Lett. 



18 



84, 48004 (2008). 

[19] J. Shao, S. V. Buldyrev, L. A. Braunstein, S. Havlin, and H. E. Stanley, Phys. Rev. E 80, 

036105 (2009). 
[20] M. E. J. Newman, Phys. Rev. E 66, 016128 (2002). 
[21] B. Bollobas, Random Graphs (Academic, London, 1985). 

[22] Fractals and Disordered Systems, Edited by A. Bunde and S. Havlin (Springer, 1996). 
[23] E. A. Leicht and R. M. D'Souza, arXiv: cond-mat/0907.0894. 

[24] L. K. Gallos, R. Cohen, P. Argyrakis, A. Bunde, and S. Havlin, Phys. Rev. Lett 94, 188701 
(2005). 

[25] A. A. Moreira, J. S. Andrade, Jr., H. J. Herrmann, and J. O. Indekeu, Phys. Rev. Lett. 102, 
018701 (2009). 

[26] H. Hooyberghs, B. V. Schaeybroeck, A. A. Moreira, J. S. Andrade, Jr., H. J. Herrmann, and 

J. O. Indekeu, Phys. Rev. E 81, 011102 (2010). 
[27] X. Huang, J. Gao, S. V. Buldyrev, S. Havlin, and H. E. Stanley, arXiv: 1010.2160. 



19 



